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We present a general numerical approach to solve the equations for bubble wall profiles in models 
with more than one scalar field and CP violating phases. We discuss the algorithm as well as several 
problems associated with it and show some profiles for demonstration found with our method. 
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For the emergence of a baryon asymmetry of the Uni- 
verse the Sakharov conditions necessarily demand devia- 
tion from thermodynamical equilibrium. This condition 
is fulfilled in first order phase transitions. They take 
place via nucleation of bubbles separating the symmet- 
ric from the broken phase. A first order phase transi- 
tion might occur at temperatures around the electroweak 
scale. It turned out that in the Standard Model (SM) 
there is no phase transition at all for Higgs masses larger 
than 72 GeV Q . Baryon number generation at the elec- 
troweak scale therefore requires more complicated mod- 
els with additional light scalar fields such as Two-Higgs- 
doublet models (2HDM), MSSM, NMSSM or further ex- 
tensions of the SM. In the MSSM there is a window for 
electroweak baryogenesis and an upper bound for the 
Higgs mass of about m# < 105 GeV with a light stop of 
mass m~ tR < m top fUJ. In the NMSSM the bounds on 
the Higgs mass are even weaker |t],|8| . 

Having established the existence of a first order phase 
transition one can start the actual calculation of the 
baryon asymmetry itself. There are several mechanisms 
described in the literature [jj-12|. All of them need 
the knowledge of the profile of the bubble wall during 
the phase transition. The kink ansatz in many situa- 
tions is a good approximation but it might be interest- 
ing to have a more refined description and to determine 
which deviations occur in the presence of potentials de- 
pending on two or more Higgs fields and one or more 
CP violating phases 13-|l6|. In fact it turns out that 
the most important value in the MSSM is the deviation 
A/3 = max v [v((3(v ) — f3(v c ))]/v c from the straight line be- 
tween the minima since the baryon asymmetry is [ pT|Jl8| 
proportional to 
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tan/3 = v 2 (x) = vf(x) + v 2 (x). Having the exact 
profile one can calculate the baryon asymmetry like flT]] 



and inve stig ate the dynamics of expanding bubbles as in 
refs. ||^||. 

To determine the bubble wall profile beyond a sim- 
ple ansatz we have to solve the equations of motion 
numerically. In the case of more than one scalar field 
this is a highly nontrivial task since simple methods like 
overshooting-undershooting which work fine in the case 
of one scalar field like overshooting-undershooting fail. 
So one has to use methods which, beginning with an 
ansatz, converge to the actual solution. They are some- 
times called "relaxation methods". Their practical im- 
plementation often is quite nontrivial. The issue of this 
letter is to present a working algorithm for the computa- 
tion of bubble wall profiles. 

We first have to find the equations of motion. In field 
theory they can be derived via Euler Lagrange equations 
from the Lagrangian density which has the general form 



C=(D^ i )+(D»<Z> i )+V($ i ,T) 



(2) 



for several Higgs fields $i (plus phases). Here D M is a co- 
variant derivative and V denotes the effective potential. 
The equations of motion also can be derived thermody- 
namically similar to Jl9| . 

Let us consider the MSSM where we may have two 
dynamical Higgs scalars (or others as a light right-handed 
stop J2|-||]) plus one CP violating phase 9. Using the real 
neutral field components 4>i, 4>2 and a relative phase 9 
the corresponding classical (tree level) Higgs potential is 



^mf 0i + ^ml4>l + mi 2 0i02 cos ( 
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In the MSSM the phase 9 also enters the stop mass 
matrix. 

To describe the phase transition we use the resummed 
one loop finite temperature effective potential (see e.g. 
ref. M) 



V T = V tree + Vi (T = 0) + Vi (T^0) + V Da 
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For the bubble wall in the MSSM without CP viola- 
tion there has been quite an interesting first numerical 
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approach to solve the problem of critical bubbles with 
two Higgs fields jL7|. In jlB) the CP profile has been in- 
vestigated in the background of a fixed Higgs profile. In 
|l6f there are detailed investigations on CP-phases with 
restriction to a straight line between the minima. 

The Euler Lagrange equations of (|j) lead to the fol- 
lowing set of coupled second order nonlinear differential 
equations 



E 1 := d^fa + Md^WO) ~ 



dVr(<f>i,h,0) 



E 2 ■= d^fo 



dv T (ct>i,<f>2,e) 
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The usual method for the SM case with only one Higgs 
field is to solve the corresponding single equation numer- 
ically by "turning around" the effective potential and 
dealing x for a time t. The problem then can be re- 
garded as an initial value problem for the negative po- 
tential and one can use an overshooting-undcrshooting 
procedure. This works well since there is only one direc- 
tion in field space. 

This situation is completely different once there are 
additional directions in field space. Again one can con- 
sider the analogous mechanical problem with the turned 
around potential. Assuming that we do not have fric- 
tion, the initial value problem is the same as trying to 
shoot a marble from one top of a hill (first minimum) 
to exactly the top of the other hill (second minimum) 
somewhere along the ridge in such a way that it comes 
to rest on the top of the second hill. Small changes in 
the initial conditions lead to a completely different shape 
of the solution. It is, in general, not possible to know 
the initial conditions with sufficient accuracy to find the 
desired solution. 

Hence we have to devise another method. We here use 
the method of minimization of the functional of squared 
equations of motion. Constraining eqs. (|))-(0) to a sta- 
tionary wall (domain wall) with velocity v w at late time 
t where the wall is already almost flat we are left with 
only one spatial dimension t perpendicular to 

the wall. Then, solving eqs. (||)-(Q) means finding field 
configurations for which 



S* = 



dx (Ef(x) + El(x) + El(x)) 



(8) 



is zero which is achieved by minimizing S3. This method 
has also successfully been used in (jlTj for the critical bub- 
ble. The approach of [p3| in principle also is a minimiza- 
tion procedure and therefore the following considerations 
are also applicable. 

Using the minimization method we have to solve a 
boundary value problem. Thus we have to use an ansatz 
for every function for which we want to find the time 
development which fulfills the boundary conditions. 
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FIG. 1. Solution (solid) for the MSSM compared to 
kink-ansatz (dashed) for tua = 450 GeV, mg = 350 GeV, 
tan/3 = 2.0. 



Our procedure works in two steps. The first step is 
to find an ansatz which is as close as possible to the 
exact solution. We will see that this is a crucial point. 
Throughout the literature we find the kink ansatz for the 
Higgs fields being quite appropriate. We therefore use for 
the N Higgs fields 
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as our ansatz as well. The Lj are determined by minimiz- 
ing (^]) . This first step is cheap in terms of computer time 
since it is only a ./V-dimensional minimization and can be 
performed very fast. Depending on the actual potential 
one may also introduce some further parameters like off- 
sets £i between the fields: cj>i — (l + tanh(-j£- + Xi) \ . 

Minimizing with respect to a few parameters is a very 
successfull first step since it reduces the actual value of 
(||) already sigificantly compared to a general function 
which only fullfills the boundary conditions. 

In figure [l] we show the solution in the MSSM for the 
simpler case where we have no CP violation (0(x) = 0) 
and the ansatz after this first step of our procedure is 
just the kink function (f>i(x) = ^-(1 + tanh(-^-)) with Li 
now determined. We can recognize that in this case the 
minimizing kink function (solid lines) is already quite 
close to the actual solution (dashed lines). The shape 
of the tunneling valley strongly depends on the CP-odd 
Higgs mass rriA ■ Small tua give a larger mixing with mn 
and a sharper bending curve. 

Unfortunately the kink ansatz does not serve as a gen- 
eral recipe. For example CP phases in general cannot be 
described very well by this ansatz and one would need 
further knowledge of their shape in the wall. For the crit- 
ical bubble with small radius or potentials with sharper 
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curves the the ansatz |9| is not appropriate any more. One 
has to choose other types of functions as ansatz and, in 
addition, the second step becomes more important. 

With this preoptimized ansatz we start the second step 
of our solving procedure, the high dimensional minimiza- 
tion. We use different numerical representations of the 
action (^) which in principle are discretizations of the 
variables over a grid. We afterwards compare the results 
obtained from different representations. In particular we 
have to discretize: 

i) The space variable x = z — v w t, 

ii) the fields 4>i(xj), 

iii) the derivatives of the fields (first and second deriva- 
tives), and 

iv) derivatives of explicitely given functions like the 
effective potential. 

ad i) We use a grid of M points Xj, j = 1 . . . M. Ac- 
cording to the special type of solution it is sometimes use- 
ful to transform the integration interval. E.g. for the kink 
ansatz which must be integrated over an infinite region it 
is useful to transform to 0i(r) with r = ^(1 + tanh(^)) 
where r e [0, 1] which is a finite interval and has the 
advantage of many interpolation points where they are 
needed, in the transition range of the kink. We have com- 
pared this to a equidistant segmentation and obtained no 
important difference in the results. By adjusting the in- 
tegration interval to the wall thickness L we take care 
that there are always enough interpolation points in the 
wall independent of L. The integration interval ranges 
between 5 and 10 times L. The regions outside of this 
interval do not contribute to (||) any more. 

ad ii) The first and quite simple method is just to store 
the function values <pi{xj) into an array. To guarantee 
smoothness here we need many grid points. We have to 
store N x M variables for all fields. 

The second approach is to interpolate the ansatz func- 
tions by smooth functions like splines. There are several 
types of splines (see refs. p4| , |25| ]). We use cubic splines, 
which are cubic polynomials through a given set of points 
with continous first and second derivative along the whole 
interpolation range. Moreover one can assign the deriva- 
tives at the endpoints of the curve which is useful for 
setting boundary conditions. For an interpolation here 
only relatively few points Xj are needed. Smoothness 
is guaranteed by definition. The price is an increase in 
computation time for each minimization step. "Nonlocal- 
ity" also has another disadvantage: If a routine changes 
one point of the curve for finding a minimum it always 
changes a larger region of the curve which may cause 
eventually larger changes in (||) than useful and maybe, 
in the worst case, even a jump into another minimum. 
We use both approaches and compare the results. 

ad iii) There are a lot of ways to discretize derivatives 
as finite difference equations connecting two or more grid 
points. We mostly use three- or four-point derivatives. 
Derivatives including more points may even impair the 
result since they require more additions and subtractions 
which arc numerically problematic. Due to quite differ- 



ent orders of magnitude in the involved numbers they 
accummulate numerical errors. For this problem see also 
[^4j,^6|. For derivative formulae see refs. p^ , |25| 

ad iv) Finally the derivatives of explicitely given func- 
tions can be calculated explicitely as well as numerically. 
This depends on the problem. For potentials like the 
MSSM potential we found it to be sufficent to differenci- 
ate numerically. 

The minimization itself is accomplished using Powell's 
method |^4| to avoid further derivatives. Minimization 
parameters are the values 4>i(xj) and 0(xj) respectively. 
So with N fields defined at M grid points we have a 
N x M-dimensional function to be minimized. Powell's 
method converges quadratically and can be used for very 
high dimensional minimizations (several hundred dimen- 
sions) . Since Powell's method is comparably slow finding 
a minimum in an almost flat valley, it might be useful in 
such cases to combine it with a Newton's method (con- 
verging quadratically as well) or the like which has a good 
convergence behaviour as well. Roughly speaking, Pow- 
ell's method finds valleys very fast and Newton's method 
finds the bottom of the minimum very fast. But the con- 
verse is not true in general. We also used downhill sim- 
plex minimization method for comparison which turned 
out to be much slowlier p4[ . Other promising algorithms 
are "leap frog" and related methods [|24|j27f| . 

Nevertheless there are still several undesired minima 
which can be categorized as follows: 

a) Some are real solutions to the equations of motion. 
As an trivial example consider 9{x) — 0, Vx which always 
is a solution to @-(0)- 

b) Fake minima due to the numerical representation of 
the functional. This is a common problem when includ- 
ing derivatives, discretized as finite differences. 

c) Minimization of (|^) is achieved by solving SS = 0. 
For a squared form S = a 2 in addition to the (desired) 
solution a = also pseudo-solutions according to 8a = 
might exist. 

There can also be combined effects between a), b) 
and c), e.g. one can have a pseudo-solution which is 
only found because of numerical "fluctuations" during 
the minimization procedure. The results often show os- 
cillations at the outer integration regions (see figure ||). 
We found that this is a combination of a) and b) as we 
see from the analytical solution of the third equation of 
motion with small 9. Then the third equation decouples 
and at small x we have the reduced equation 

6" + -6' - \m\ 2 1 tan (36 = (10) 

which has the solution 

6(x) = exp{-2x/L} cos(lux/L), (11) 

a damped oscillation with frequency 
u> = y4 — L 2 \m 2 2 \ tan/3 and and growing amplitudes for 
x — > —oo. 
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FIG. 2. Oscillation of an undesired solution of 6(x) near 
the origin,where 6 is small and decouples from <j>\ and <f>2- 

This is what we find numerically. But since our bound- 
ary condition were 0'(±oo)=O and we took as an ansatz 
the zero function 0(x) = 0, we find this behaviour to 
be caused by the minimization routine. The routine 
came through a configuration around 8(x) = 0. It dif- 
fers inessentially from the ansatz but together with the 
chosen boundary conditions this configuration fits into 
the oscillation described above. Without this pseudo- 
solution 6(x) would remain zero for all X ; CIS intended. 
Fortunately in the region of interest within the bubble 
9(x) remains zero since the chosen parameters permit no 
CP violation there. This is typical for situations where 
an angle with vanishing moduli is not well defined any 
more. 

Altogether this implies the possible problem that, 
starting from an ansatz which is in the vicinity of such an 
apparent solution, the algorithm might never converge to 
the desired exact solution. Therefore the importance of 
a good ansatz is obvious. 

What can be done to avoid some of the problems and 
rate the quality of a minimum found? First we see that 
it is important to have a well prepared ansatz as near as 
possible to the desired solution to avoid reaching an un- 
wanted local minimum after the time consuming second 
minimization step. As long as we have energy conserva- 
tion, as in our example, one can check the quality of the 
results easily. It must be T — V = or T/V = 1 and 
one can check the deviation. Here V is the QFT effective 
potential, the mechanical analog has the opposite sign, T 
is the kinetic energy density. It is possible to get results 
where T/V deviates only by few percent from 1.0. The 
energy check is already quite appropriate, it strongly de- 
pends on the quality of the solutions. It is not useful to 
add the energy conservation, which would be fullfilled by 
the real solution automatically, since it is another min- 
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FIG. 3. Unphysically small values of rriA ~ O(10 GeV) 
demonstrate the deviation of the solution (solid) from the 
straight line (dashed). 

imum finding task with all problems we have described. 
Only in the stationary case T — V could actually be used 
instead of (||) whereas our method still holds in a more 
general case. One can think of further identities that have 
to be fulfilled, as in [p3| . The described special behaviour 
of unwanted oscillations can be suppressed by restricting 
the parameter space of the minimization procedure at the 
boundaries. One can also increase the number of steps in 
the procedure to improve the precision step by step for 
the cost of computing time. 

Since the kink solution is a good approximation to the 
bubble wall profile equations in a lot of models, we com- 
pared our result with |l7j . In that paper the radial sym- 
metric equations for the critical bubble were solved. Thus 
the solutions have a different shape but in field space 
there should be the same behaviour concerning the de- 
viation from the straight line between the minima. And 
we indeed can confirm the results. A/3 is in the range 
of 1CP 3 to 1CP 2 for rriA around several hundred GeV 
H[l7|]. Only for (experimentally excluded) small values 
of rriA ~ O(10GeV) we find a considerable deviation from 
a straight line demonstrating a major deviation from the 
ansatz due to the minimization routine (see fig. [5]). 

CONCLUSIONS AND PROSPECT 

In order to solve the phase transition problem for find- 
ing bubble wall solutions we presented a working algo- 
rithm. Often a simple ansatz like the kink turns out to 
be a good approximation. Especially in the 2HDM and 
the MSSM the deviations from this ansatz are small at 
least in some expressions like the surface tension. But for 
example for the baryon number asymmetry or the CP 
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violating phase this is not sufficient since these indeed 
depend on the small deviation from the simple ansatz: 
The latter would give zero for A/3. In some models this 
would prevent baryogenesis at all. 

With our method we can confirm the results of JlTj] 
for the one loop MSSM effective action. The methods 
presented here can also be extended to investigate the 
dynamics of bubble expansion in more detail. For this 
we have to go beyond a constant critical temperature and 
include reactions of the bubble wall with a background 
fluid of particles as in |21|. These questions have to be 
studied further Q. An investigation of the NMSSM 
with three Higgs fields and a stronger deviation from a 
straight line in field space is in progress as well. 
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